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This work measures the time to equihbrium for the multicanonical method on the 2D- 
■ ■ Ising system by using a new criterion, proposed here, to find the time to equilibrium, 

^ ' teq, of any samphng procedure based on a Markov process. Our new procedure gives the 

O ' same results that the usual one, based on the magnetization, for the canonical Metropolis 

O , sampling on a 2D-Ising model at several temperatures. For the multicanonical method 

we found a power-law relationship with the system size, L, of teq=0.27(15)-L^ *'''^^' 
and with the number of energy levels to explore, ^e, of teq=0.7(13)fcg , in perfect 
agreement with the result just above. In addition, a kind of critical slowing down was 
observed around the critical energy. Our new procedure is completely general, and can 
' be applied to any sampling method based on Markov processes. 

Keywords: Multicanonical method, time to equilibrium, Markov processes, flat- 
04 ' histogram algorithms 
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. 1. Introduction 

P 

The multicanonical method, first proposed by B. Berg and T. Neuhaus ^'^ and 
^ also known as entropic sampling is a powerful method to explore the phase space 

of many physical systems on broad energy ranges. It has been sucessfully employed, 
' not only to the study of both first and second order phase transitions in statistical 

^1 , mechanics but also to a plenty of many other areas, from optimization problems 

in computational sciences ^ to protein folding ^'^'^o and even covering rules for 
quasicrystals (for a review on the applications of the multicanonical method, see 
^^). This method uses a Markov process with an extended Metropolis accepting 
ratio to sample any state x in the phase space with a probability P{x) proportional 
to the inverse of the density of sates g{E), 



(1) 
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with E{x) the energy of state x. Then, the Metropolis accepting ratio for a move 
from X to x' becomes 



By this method, the histogram of samples on the energy axis, V{E), is flat, so long 
g{E) is the actual density of states. If not, a better approximation gi{E) can be 
found by dividing V{E) into P{x), i.e., 



The multicanonical method is actually an iterative procedure that, first, proposes 
some trial function g{E) and, then, improves it by successive applications of the 
procedure described above. 

As every Markov process does, the multicanonical method needs a number tgq 
of steps before reaching the desired probability distribution, that is, before starting 
to take the samples. For some sampling distributions, like the Boltzmann one, it is 
easy to find narrowed- distributed quantities, like the magnetization, such that teq 
can be measured just by tracking the process until that quantity oscillates around 
a mean value and inside a 10% confidence interval, for instance. That is what is 
called statistical control. The well-established procedure of computing the non- linear 
correlation function and the non-linear correlation time is just a sophistication 
of this procedure, which gives good results if the quantity evolves to its equilibrium 
distribution as an (exponential decay. For the multicanonical method, in contrast, it 
is not so easy to imagine a narrow-peaked quantity that can be used for this purpose. 
The typical procedure to find teq is to run two or more instances of the system, and 
to keep track in all of them of some quantity of interest, i.e. the magnetization, 
until all instances join together. The time to equilibrium teq is taken as the number 
of steps, measured in Monte Carlo steps per site (MCSS), before this joint. Despite 
the fact that this second criterion is easy to improve, it deserves the statistical 
precision of the first one. 

This paper introduces a new criterion to measure t^q and uses it to investigate 
the time to equilibrium of the multicanonical sampling method. The new criterion 
is based on the chi-squared (x^) deviation between the desired limit distribution 
and the actual energy distribution after t steps. First, we introduce this chi-squared 
procedure, and we test it against the more traditional one, based on the magneti- 
zation, for the usual Metropolis sampling with Boltzmann distribution. Second, we 
use the chi-squared procedure to investigate the time to equilibrium for the multi- 
canonical method on 2D-Ising systems, and how it depends on the system size and 
the energy. Finally, we summarize our main results and conclusions. 

2. The Procedure 



A{x'\x) = mm{l,g[E{x)]/g[E{x')]} 



(2) 



g^{x) c< V{E)g{E) 



(3) 



Consider a sampling method based on a Markov process that takes any state 
X with some probability P{x), which is function of E{x) alone. Starting from 
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some initial distribution Po{x), the Markov process transforms it at every step 
until it reaches the desired distribution P{x). Let us run n times the Markov 
process, starting form the same initial distribution Po{x) but running with different 
seeds of the random number generator. Look at the energy the system has after 
t steps on each run, and cumulate the histogram Vt{E) from these energy values. 
If Vt{E) actually corresponds to the limit distribution P{x), it should approximate 
the expected probability distribution for E, P{E)=g{E)P{x). The deviation is 
a good statistics to verify this point. It is given by 

where the sum extends over all energy values with more than five or ten samples. 
For large n values, obeys a chi-square distribution of r=kE — 1 degrees of freedom, 
with probability density 

2 (x2)(''-2)/2e-xV2 
-^(^ ) = r(r/2)2'-/2 ' 

where fc^; is the number of energy levels in the sum (^^ . This distribution has 
population mean /x=r and standard deviation a=y/2r. If we repeat m times the 
whole study and look at the average x^ of that m values, x^ distributes like a 
Gaussian, with fj,=r and a=i/2r/m. That is, 95,45% of the values of x^ lies between 
r — 2y/2r/m and r — 2-\/2r/m. This is only true if Vt{E) corresponds to the limit 
distribution P{E). That is what is called statistical control 

So, the new procedure to find the time to equilibrium is just to compute x^; as 
described above, and to look at the time, teq, when the statistical control starts. 
Figure 1 shows how this criterion works for a Metropolis sampling with limit Boltz- 
mann distribution on a 8x8-Ising model. The temperature was r=3.0, and the 
system always started from the state of minimal energy. In this example, only 
kE=22 energy levels were taken into account to compute x^(i) (Eq- (4)), and the 
whole study was done m=10 times to compute x^) with n=10'* runs each. There- 
fore, statistical control is reached when x^ oscillates around r=21, and between 
r - 2^/2r/■m=^6.9 and r + v^2r/m=25.1. We found teq = 33(3). 

Figure 2 shows t^q for the same 8x8-Ising model, when sampled by the Metropo- 
lis method with limit Boltzmann distribution at several temperatures. Two criteria 
were used to measure teq- the x^-criterion described above and another one based 
on the statistical control on the magnetization M, as follows: we looked at the 
magnetization the system had after t steps on each one of the n runs and computed 
their mean value <M{t)>. The time to equilibrium was chosen as the time when 
statistical control started, with ii—<M>t (the average magnetization at tempera- 
ture T) and a—y^Txs{T)/m (with Xs{T) the magnetic susceptibility at temperature 
T, Xs{T)=l/kBT[<M^>T - <M>^]). All other parameters in the simulation were 
the same as in figure 1. We can see an excellent agreement between both criteria 
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Fig. 1. Example of the x^-criterion to compute the time to equilibrium teq for a Metropolis 
sampling method with limit Boltzmann distribution (T=3.0) on a 8x8-Ising model. 



for all temperatures. That is, our new criterion works well by measuring the time 
to equilibrium of this Markov process. 




Temperature 

Fig. 2. Time to equilibrium teq for a Metropolis sampling method with limit Boltzmann distribu- 
tion at several temperatures on a 8 x 8-Ising model, computed both with the new x^-criterion and 
by looking at the statistical control on the magnetization. The two criteria agree for all sampled 
temperatures. 



3. The multicanonical case 

With the new x^-critcrion on hand, wc investigate the time to equilibrium for 
the multicanonical sampling method on two-dimensional Ising models of several 
sizes. For all simulations the actual density of states g{E), computed by P.D. Beale 

was used from start to perform every multicanonical step. In all cases we chose 
n=10* and m=10, as before. 

Figure 3 shows time to equilibrium against system size, when the whole energy 
range of positive temperatures —2L'^<E<0 is sampled. Through linear regression. 
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we obtain a power-law behavior, icq=0.27(15)L^'^''*^"'^^'' f with correlation coefficcients 
r=0.9988 and i^=2416.72. That is, the time to equihbrium increases with system 
size likeL2-80(i3). 




Fig. 3. Dependence of the time to equilibrium teq on system size L for the multicanonical sampling 
method on LxL-Ising models. The true density of states g{E) was used from start at every 
multicanonical step. 



Next, we look at how the time to equilibrium changes with the number of en- 
ergy levels (the number of bins in the histogram) sampled by the multicanonical 
method. For that purpose a set of multicanonical simulations was performed on 
Ising models of system sizes L—6, 8, 10, 12 and 16. All energy intervals start from 
Einf=—^L'^ and extend up through ks energy levels. Figure 4 shows teq against ks 
for all system sizes. We obtain, again, a power-law behavior, feg=0. 70(13) 
with correlation coefficients r=0.9783 and i^=598.19. Summarizing, tgq grows like 

,1.40(11) 

Finally, we hold constant the relative length of the energy interval and we inves- 
tigate how the time to equilibrium changes with the position of this interval along 
the energy axis. With Emin=^'^L'^ , an energy interval of length 0.16|i?TOm| was set 
at different places on the energy axis, actually around approximate energy values of 
0.92, 0.84, 0.75, 0.68, 0.59, 0.52, 0.43, 0.36, 0.25 and 0.10 times Emm- The system 
reached the energy interval from below by starting at a configuration of minimal 
energy and by accepting every random flip that increases the energy. This proce- 
dure fix the initial distribution Po{x). A whole set of simulations was performed 
on systems of sizes i=6, 8, 10, 12 and 16, as before. Figure 5 summarizes these 



results. We try a finite-size scaling around the critical energy ii^c=0.7540328£'„iCT 
(the mean energy at the critical temperature Tc~2. 269185 for an infinite-size sys- 
tem) by plotting L^^-^Heq against LP '^{E — Ec)/Ec. Although the scaling is not 
"All error bars in this paper are two-sigma. 
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Fig. 5. Time to equilibrium teq for the multicanonical samplig method on energy intervals of length 
0.16|Bmin| (with E^j„=— 2L^), placed at different places on the energy axis. The horizontal axis 
represents the value for the central energy E in the interval. A certain kind of critical slowing down 
can be observed around a critical energy of Bc=0.7540328B^i„. The simulations were performed 
on two-dimensional Ising systems of sizes L=6, 8, 10, 12 and 16. 



a perfect one, we observe a certain kind of critical slowing down, i.e. the time to 
equilibrium teq grows if the energy interval lies near the critical energy and diverges 
with increasing system size. 

4. Conclusions 

As every Markov process does, multicanonical sampling needs some steps be- 
fore reaching the desired sampling distribution. We have found that this time to 
equilibrium (measured in MCSS) changes mainly with the number of energy bins, 
kE, to be sampled. Our main result is that this dependence obeys a power law, 
feg=0. 70(13) A:^''°'^^\ at least for the two-dimensional Ising model. If one wants 
to sample the whole energy interval of positive temperature for several system 
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sizes (which consists of kE=L'^/'2 energy bins), this resuh predicts that teq will 
grow with system size, L, as teq=0.27(5) i^'^"*^^^-*. A direct measurement gave us 
feq=0.27(15)L-^-^*'(^^\ in perfect agreement with the prediction above. 

This result has important consequences. First, it supports the usual procedure 
with multicanonical samplings of dividing the energy interval of interest in many 
small intervals and sampling each one separately. This is not just to parallelize the 
method, it actually reduces the total time to equilibrium, because it grows faster 
than linear with ks- Second, it says us that we have to be much more careful by 
performing multicanonical samplings on more than one parameter. If we construct 
an histogram on two parameters (let us say, the energy and the magnetization, for 
instance) with k bins on each one, the time to equilibrium will grow like and 
that can be much faster than expected (so far our results can be extended to this - 
somewhat different - situation). 

By dividing the energy range of interest in a fixed number of small intervals and 
by using a multicanonical sampling on each one, we have found that the equilibration 
times is larger for energy intervals around the critical energy (the mean energy 
at critical temperature), and that those times grow with system size like L'^-°, 
approximately. Although the finite-size scaling is not perfect, it tells us that we 
have to be more careful on energy intervals around the critical energy, and that the 
time to equilibrium should be measured just there, if we want to assign the same 
number of MCSS equilibration steps for all intervals. 

All these results were obtained by using a new criterion to estimate the time to 
equilibrium of any sampling Markov process. We have shown that this new criterion 
gives the same equilibration times than the usual one (based on the magnetization) 
for the canonical Metropolis sampling of a 2D-Ising model at several temperatures. 
Although general, this new method is especially useful for sampling procedures that 
construct a flat histogram on the energy axis, like the multicanonical method, the 
flat-histogram method of J.-S. Wang the sampling method of F. Wang and D. P. 
Landau or the broad-histogram sampling of P.M.C. de Oliveira We are sure 
that our criterion will work properly for all those sampling procedures, revealing 
important issues of its behavior that can encourage its proper use and improve its 
performance. 
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